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\  ABSTRACT 

A  statistical  model  is  analyzed  for  the  growth  of  bubbles  in  a 
Rayleigh-Taylor  unstable  interface.  The  model  is  compared  to  solu¬ 
tions  of  the  full  Euler  equations  for  compressible  two  phase  flow, 
using  numerical  solutions  based  on  the  method  of  front  tracking.  The 
front  tracking  method  has  the  distinguishing  feature  of  being  a 
predominantly  Eulerian  method  in  which  sharp  interfaces  are 
preserved  with  zero  numerical  diffusion.  Various  regimes  in  the  sta¬ 
tistical  model  exhibiting  qualitatively  distinct  behavior  are  explored.  It 
appears  that  the  parameters  in  the  statistical  model  can  be  set  from  first 
principles  on  the  basis  of  comparison  with  numerical  solutions  of  the 
full  Euler  equation.  ^ 

Part  of  this  work  was  performed  during  a  visit  at  the  Cornel]  Center  for  Theory  and  Simulation 
Permanent  Address 
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A.  Results,  Methods  and  Goals 

The  purpose  of  this  paper  is  to  formulate  a  program  for  the  statistical  analysis  of 
the  long  time  growth  of  bubbles  formed  by  Rayleigh-Taylor  unstable  interfaces.  In 
addition,  this  paper  will  carry  out  some  initial  stages  of  the  program.  There  are  two 
essential  ingredients  to  this  program.  The  first  is  an  unpublished  statistical  model1 
for  bubble  growth  and  cannibalization,  developed  by  one  of  the  authors  (D.  H.  S.) 
and  J.  A.  Wheeler  in  1961.  An  important  prediction  of  this  model  which  we  further 
investigate  and  quantify  is  that  bubble  merger  leads  to  an  accelerated  growth  rate  of 
the  interface.  The  second  ingredient  is  a  front  tracking  hydrodynamics  code,  suitable 
for  the  study  of  unstable  interfaces,  developed  by  the  authors  and  coworkers.2-7 

There  are  two  main  new  results  in  this  paper.  We  have  carried  out  a  careful 
study  up  to  late  times  of  the  dynamics  of  a  single  mode  (bubble  and  spike)  for  a 
Rayleigh-Taylor  interface  in  the  compressible  regime.  Moreover  we  have  identified 
new  regimes  and  correlations  in  the  interactions  of  ensembles  of  large  numbers  of 
bubbles  through  a  study  of  the  statistical  model  for  bubble  dynamics.  These  results 
are  placed  in  a  context,  which  if  developed  further,  could  allow  a  first  principles 
understanding  of  the  Rayleigh-Taylor  hvixing  layer. 

The  statistical  model  contains  two  free  parameters,  which  are  the  velocity  of  a 
single  bubble  and  the  relative  location  for  merger  of  two  adjacent  bubbles.  These 
parameters  can  be  viewed  as  describing  the  essential  features  of  the  one  and  two 
body  bubble  dynamics.  By  a  study  of  one  and  two  body  bubble  dynamics  in  the  full 
Euler  equations  of  fluid  dynamics,  using  the  front  tracking  code,  we  can  set  these 
parameters  and  analyze  their  dependence  on  density  ratios,  compressibility,  equation 
of  state,  geometry,  viscosity  and  surface  tension.  In  principle  their  dependence  on 
spatial  dimensions  can  similarly  be  determined,  by  use  of  three  dimensional  compu¬ 
tations. 

Once  the  parameters  of  the  statistical  model  have  been  fixed,  it  can  be  studied 
on  its  own  right.  The  most  basic  question  is  the  possible  existence  of  distinct  parame¬ 
ter  regimes  or  domains  of  initial  data  producing  qualitatively  distinct  behavior  (e.g. 
regular  or  quasi-uniform  vs.  runaway  or  single  finger  dominant  modes).  Any  well 
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defined  boundaries  or  separatrices  which  might  exist  between  such  domains  could  be 
thought  of  as  defining  a  phase  transition  in  the  statistical  behavior  of  an  ensemble  of 
bubbles.  Within  a  single  solution  domain,  the  time  evolution  of  statistical  averages 
such  as  the  average  bubble  size,  height  or  velocity  are  important.  Moreover  the  sta¬ 
bility  or  insensitivity  of  these  quantities  to  statistical  noise  in  the  initial  conditions 
will  also  be  important. 

The  program  presented  here,  for  the  analysis  of  bubble  dynamics,  is  only  part  of 
a  larger  question,  which  is  a  theory  of  the  Rayleigh-Taylor  mixing  layer.  The  mixing 
layer  can  be  analyzed  in  terms  of  various  regions,  of  which  the  bubble  regime  is  only 
one,  located  at  the  edge  of  the  mixing  layer  closest  to  the  undisturbed  heavy  fluid. 
There  is  a  second  regime,  dominated  by  spikes  or  droplets  at  the  edge  of  the  undis¬ 
turbed  light  fluid.  A  theory  known  as  dusty  gas  could  be  appropriate  for  this  region. 
It  has  not  been  determined  whether  additional  interior  regimes  are  needed  for  a  full 
description  of  the  mixing  layer.  Such  a  theory  should  describe  blobs  of  heavy 
material  moving  at  a  terminal  velocity  and  shedding  vortices  (Atwood  ratio  A  <  1)  or 
in  free  fall  without  vortices  (A  =  1),  while  the  pinch  off  of  sheets  into  spikes  and 
spikes  into  droplets  would  seem  to  depend  on  surface  tension.  Regimes  in  the  time 
evolution  of  the  mixing  layer  can  also  be  identified,  such  as  the  small  amplitude 
(linear),  the  coherent  individual  mode,  the  mode  competition  and  the  chaotic 
regimes.8  Mixing  theories  have  been  advanced  in  a  number  of  other  contexts  as 
well  9-11  See  ref.12  for  a  previous  study  of  the  compressible  Rayleigh-Taylor  prob¬ 
lem. 

B.  Bubble  Merger  and  Acceleration  of  the  Bubble  Interface 

The  present  study  goes  beyond  previous  Rayleigh-Taylor  studies  in  a  number  of 
respects.  In  doing  this,  we  make  use  of  the  unique  capabilities  of  the  front  tracking 
code.  The  detailed  study  of  the  one  body  bubble  problem  presented  here  for 
compressible  Rayleigh-Taylor  fingering  seems  to  be  the  first  attempt  to  set  the 
parameters  of  the  Sharp-Wheeler  model  from  a  first  principles  calculation.  Youngs 
has  argued8  that  deeply  nonlinear  multi-mode  computations  of  the  chaotic  regime  are 
necessary  to  understand  the  Rayleigh-Taylor  chaotic  mixing  regime.  Since  such  com¬ 
putations  are  at  best  only  marginally  possible  and  of  uncertain  validity  in  their  predic¬ 
tion  of  statistical  mixing  phenomena,  it  is  significant  that  some  statistical  information 
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may  be  extracted  from  the  study  of  the  one  and  two  bubble  problems.  A  test  of  these 
predictions  would  be  to  compare  them  to  direct  simulation  of  the  deeply  nonlinear 
chaotic  regime,  using  front  tracking  and  other  methods.  The  computational  analysis 
of  the  Sharp-Wheeler  model  is  also  new.  The  identification  of  a  transient  exponen¬ 
tially  growing  runaway  mode  in  the  bubble  merger  process  and  the  importance  of 
nearest  neighbor  correlations  in  limiting  bubble  amalgamation  appears  to  be  new. 
The  constant  acceleration  rate1,8,13  of  the  bubble  interface  has  been  confirmed  in  the 
examples  of  initial  data  considered,  but  the  rate  has  been  found  to  depend  on  initial 
conditions,  and  seems  not  be  a  universal  constant  at  least  for  the  time  periods 
explored  in  this  paper.  Youngs’  argument8  in  favor  of  universality  was  based  on  a 
scaling  argument  and  on  an  absence  of  extra  length  scales  in  the  problem.  Our  exam¬ 
ple  of  non-universal  acceleration  rates  depends  on  introducing  an  extra  length  scale. 
This  was  introduced  by  slightly  varying  the  bubble  radii,  so  that  the  deviation  in 
radius  was  much  smaller  than  the  mean  radius.  Because  of  the  presence  of  a  length 
scale  in  the  initial  data,  we  do  not  contradict  the  scaling  argument  of  Youngs,  but  we 
do  suggest  that  his  hypothesis  of  the  absence  of  extra  length  scales  may  be  valid  only 
for  a  restricted  class  of  Cauchy  data  or  for  extremely  late  times.  Experimental 
results  on  carefully  prepared  interfaces  reveal  a  universal  behavior,  with  the  bubble 
interface  height  h  =  h(t)  satisfying  the  relation  h  -  agAt2  and  a  ~  .06  (two  dimen¬ 
sion)  or  a  ~  .07  (three  dimension)  15,  where  A  is  the  Atwood  ratio,  A  =  ^  +  . 

Here  pa  is  the  density  of  the  gas  above  the  interface,  and  p *,  is  the  density  of  the  gas 
below  the  interface. 


Other  mechanisms  leading  to  a  constant  acceleration,  agAt2,  of  the  bubble  inter¬ 
face  are  possible.  For  an  initial  time  period,  the  bubbles  accelerate  (are  in  free  gravi¬ 
tational  rise).  For  the  strongly  compressible  case  considered  here,  the  bubbles  attain 
velocities  up  to  Mach  2  relative  to  the  sound  speed  of  the  heavy  gas.  Comparing  our 
values  for  free  rise  renormalized  gravity  (~  .5 Ag)  for  bubbles,  we  get  a  position 
h  ~  .25 Agt2  by  this  mechanism,  which  is  approximately  four  times  the  position 
predicted  by  the  merger  mechanism,  according  to  Read’s  data.  Free  rise  describes  an 
earlier  time  period  than  the  merger  process,  and  compressibility  plays  a  major  role  in 
setting  the  duration  of  this  earlier  regime. 


A  second  possible  mechanism  for  constant  acceleration  is  that  fandom  initial 
data  may  contain  a  mixture  of  large  and  small  wavelengths.  The  small  wavelengths 
grow  more  rapidly  at  first  but  then  saturate  and  eventually  the  initially  slower  large 
wavelengths  become  faster  and  win  out.  Thus  large  structures  may  be  latent  in  the 
initial  data,  and  emerge  gradually  over  time,  leading  to  an  acceleration  of  the  bubble 
interface.  Neither  of  these  alternate  mechanisms  appear  to  be  important  in  Read’s 
data. 

Portions  of  our  results  have  been  announced  previously.15-17  (See  ref.18  for  a 
further  discussion  of  this  problem). 

II.  A  FORMULATION  OF  THE  STATISTICAL  MODEL 

The  main  idea  in  the  statistical  model  is  to  introduce  a  very  reduced  set  of 
parameters  ("idealized  interfaces")  and  a  simplified  interface  dynamics,  which  cap¬ 
tures  the  essential  features  of  the  full  flow  in  the  mixing  layer.  For  comparison  pur¬ 
poses,  we  also  need  a  map  from  arbitrary  to  idealized  interfaces.  It  is  then  hoped 
that  the  stable  statistical  properties  (average  bubble  size,  etc.)  of  the  model  will 
correspond  to  stable  statistics  of  the  full  flow  and  that  after  adjustment  of  model 
parameters,  the  dynamics  of  such  stable  statistical  properties  will  approximate  the 
true  dynamics  of  such  stable  properties  of  the  full  flow. 

The  first  approximation  we  make  in  going  from  the  full  flow  to  the  model  is  to 
assume  that  the  interface  position  at  time  t  is  a  single  valued  function  z  =  z(x, y,r)  of 
the  point  x,  y.  We  do  not  assume  r,  y,  :  to  be  a  rectangular  set  of  coordinates,  and 
so  various  curvilinear  geometries  are  included  in  this  analysis.  Let  us  suppose  that 
the  heavy  fluid  is  located  in  the  region  z  <  z(x,y,t)  below  the  interface  and  that  the 
accelerating  force  is  in  the  direction  of  the  positive  z  axis. 

Since  a  Rayleigh-Taylor  unstable  interface  at  late  time  will  be  far  from  single 
valued  and  since  we  are  developing  here  a  theory  of  the  bubble  region  only,  within 
the  mixing  layer,  we  regard  z(x,y,t )  as  representing  a  "bubble  envelope"  equivalent 
to  the  true  interface,  defined  precisely  by  the  imposition  of  mass  conservation  along 
each  line  x  =  const.,  y  =  const.  Thus  if  the  entire  multivalued  interface  lies  below 
some  reference  height  zq,  we  define 


where  L  is  the  length  of  the  portion  of  the  line  x  =  const.,  y  =  const,  located  in  the 
light  fluid  above  zq. 

Next  we  fix  a  set  of  bubble  boundaries,  which  in  mathematical  language  is  a  til¬ 
ing  of  the  surface  r  =  0  by  polyhedra.  In  two  space  dimensions,  the  bubbles  thus 
define  a  decomposition  of  the  x  axis  into  non-overlapping  intervals.  Even  in  this 
simplest  case,  we  do  not  propose  a  general  algorithm  for  the  placement  of  bubble 
boundaries,  but  we  will  return  to  this  question  as  part  of  a  detailed  analysis  of  when 
or  whether  merger  occurs  in  the  two  body  bubble  interaction  problem. 

Finally  the  model  assumes  that  the  fluid  is  piecewise  constant,  and  constant 
within  each  bubble.  Given  a  bubble  decomposition  (polyhedral  tiling)  of  a  general 
single  valued  interface  z  =  r(.r,y,r),  the  principle  of  mass  conservation  again  defines 
a  corresponding  interface  within  the  model,  i.  e.  piecewise  constant  on  each  bubble. 

The  kinematics  of  the  model  is  now  determined;  we  must  specify  its  dynamics. 
There  are  two  processes  allowed  in  the  dynamics  of  the  model:  vertical  motion  and 
merger.  We  assume  that  each  bubble  moves  vertically  with  a  velocity  z  which  is 
independent  of  the  bubble  height  and  the  other  bubbles.  Thus  i  is  a  function  of  the 
bubble  shape  alone.  We  note  that  for  a  circular  bubble  of  vacuum,  with  radius  r,  ris¬ 
ing  in  an  incompressible  fluid  in  three  space  dimensions,  Taylor’s  formula  gives 

l 

-i=.48(gr)T.  (1) 

In  two  dimensions,  i  =  z(r),  is  thus  assumed  to  be  a  function  of  the  bubble  radius 
alone.  As  a  further  approximation  in  three  dimensions,  we  define  an  effective  radius 
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r  =  (  2^  ) 2  and  postulate  that  z  =  z(r )  in  this  case  also.  In  effect,  this  means  that 
the  bubbles  do  not  depart  too  greatly  from  a  circular  cross  section. 

The  remaining  aspect  of  the  model  bubble  dynamics  is  merger.  It  is  postulated 
that  adjacent  bubbles  merge  when  one  is  sufficiently  far  ahead  of  the  other.  After 
merger,  there  is  a  single  bubble,  with  cross  section  equal  to  the  union  of  the  cross 
sections  of  its  two  constituents  and  a  height  z  determined  by  conservation  of  mass 
within  this  cross  section.  If  the  heights  and  radii  of  the  two  bubbles  before  merger 
are  z,  and  r,-,  i  =  1,2,  then  we  formulate  the  merger  criterion  as 

Z2  -  z i  S  —  m(rj)  (2) 


lm \Ar^  *>  aw  v.  < 


where  n  s  ri  and  m(r)  is  an  unknown  function  of  r. 

Thus  the  statistical  model  is  defined  by  two  unknown  functions,  i(r)  and  m(r), 
which  express  the  essential  features  of  the  one  and  two  body  bubble  problems  respec¬ 
tively.  These  functions  must  then  be  set  either  by  some  external  theory  or  experi¬ 
ments  or  by  a  numerical  study  of  these  one  and  two  body  problems  in  the  full  Euler 
equations. 

One  aspect  of  this  model  is  that  after  merger  the  new  bubble  may  also  satisfy 
the  merger  criterion  with  an  adjacent  bubble  thus  giving  rise  to  a  cascade  of  mergers 
at  a  single  instant  of  time.  For  example,  a  cascade  occurs  if  the  next  nearest  neigh¬ 
bor  of  a  large  bubble  is  smaller  than  its  nearest  neighbor  and  the  height  of  the  next 
nearest  neighbor  is  less  than  the  height  of  the  nearest  neighbor.  Thus  the  merger 
process  in  this  model  for  some  distributions  of  bubbles  consists  of  more  than  two 
body  interactions. 

III.  A  QUALITATIVE  STUDY  OF  THE  STATISTICAL  MODEL 

We  assume  a  probabilistic  distribution  of  initial  data  in  the  form  of  bubble 
heights  and  widths.  It  should  be  noted  that  the  model  of  Section  2  is  actually  deter¬ 
ministic,  and  becomes  a  statistical  model  at  this  point  alone.  For  simplicity  we  now 

consider  bubbles  in  two  dimensions.  From  a  fundamental  point  of  view,  we  imagine 

/. 

a  translation  invariant  ensemble  of  random  bubble  interfaces,  constrained  to  be 
defined  by  single  'alued  piecewise  constant  functions  z  =  z(x).  This  ensemble  is  fully 
described  by  a  set  of  n-point  functions  or  correlation  functions,  i.e.  conditional  pro¬ 
bability  distributions  for  a  sequence  of  n  adjacent  bubble  widths  and  heights.  The 
model  dynamics  is  equivalent  to  an  infinite  system  of  coupled  partial  differential 
equations  for  these  n  point  functions.  This  system  does  not  close,  in  the  sense  that 
the  equations  for  the  j  point  functions,  lS/Sn,  depend  on  the  n  +  1  point  func¬ 
tion.  We  invoke  the  simplest  closure  hypothesis,  assuming  that  the  statistical  proper¬ 
ties  of  points  are  independent,  or  equivalently  that  the  truncated  two  point  function  is 
zero.  In  the  next  section  it  will  be  shown  that  this  closure  hypothesis  excludes  an 
important  range  of  phenomena.  Let  T|(w,r,r)  be  a  density  function  for  the  distribu¬ 
tion  of  bubble  widths  w  and  heights  z  at  time  t.  The  closure  hypothesis  leads  to  a 
nonlinear  set  of  coupled  partial  differential  equations  for  t).  Namely 


Here  B  and  D  are  the  birth  and  death  terms  respectively,  resulting  from  bubble 
merger,  while  the  left  hand  side  expresses  the  noninteracting  growth  of  a  single  bub¬ 
ble.  Specifically 


B(w,z,t )  =  £j\(w',z,,t)T\(w",z",t)  | z(w’)  -  z(w")\dw', 

D(w,z,t)  =  2rK>v,z,r)£n(w',r',f)  |i(w’)  -  z(w)\dw', 
where  T%  and  Tb  specify  domains  of  integration,  defined  as  follows: 

—  D  width  O  Dmass  O  Emerge 

rb  =  {h  ',  z':  0  s  w  s  w1  and  z'  —  z  =  m(w)} 

U  {w',  z':  0  :£  w'  ^  w  and  z  —  z'  =  m(vv')}  . 

Here  DwUith,  Dmass  and  Dmergt  denote  the  conservation  of  width  and  mass  and  the 
merger  height  constraints  respectively. 

Dwid,h  -  {'v\  z',  w",  z":  w'  +  w”  =  w  and  w',  w"  2:  0}, 

Dmass  =  {'v',  Z',  W",  z":  w'  z'  +  w"  Z  '  =  W  z}, 

=  {w',  z',  w",  z":  w"  £  w'  and  z'  -  z"  =  m(w")} 


U  {u'',  z',  m’'',  z":  w'  <  w"  and  z"  —  z'  —  m(w')} 

In  order  to  continue  the  analysis  we  assume  a  distribution  in  the  heights  z  and  bub¬ 
ble  widths  w,  which  are  characterized  by  time  dependent  parameters,  such  as  their 
means  and  covariances.  Differentia]  equations  can  then  be  derived  for  these  parame¬ 
ters,  which  however  will  not  close,  in  the  sense  that  time  derivatives  of  the  distribu¬ 
tion  will  not  in  general  lie  in  the  same  finite  dimensional  subspace  characterized  by 
these  parameters.  Again  we  impose  a  closure  hypothesis,  which  corresponds  to  a 
projection  of  the  dynamics  onto  a  fixed  finite  dimensional  subspace.  We  do  not  pur¬ 
sue  this  point  of  view  further  in  the  present  paper. 


IV.  A  NUMERICAL  STUDY  OF  THE  STATISTICAL  MODEL 

There  are  two  constants,  cj  and  c2,  which  characterize  the  statistical  model, 
namely 

1 

-i(r)  =  ci(jr)7  (3) 

and,  with  the  further  simplifying  assumption  that  m(r)  is  linear, 

m(r)  =  c2r,  (4) 


cf.  (1)  and  (2).  To  begin  with,  we  consider  the  model  in  the  Jt,  z  plane.  Then  the 
mean  bubble  size  fixes  the  x  units,  c 2  fixes  the  relation  of  the  z  units  to  the  x  units, 

c } 

and  g  fixes  the  relation  of  the  r  unit  to  the  time  units.  —  is  a  dimensionless  free 

c2 

parameter  of  the  model.  Moreover,  the  mean  height  fixes  the  origin  of  the  z  units. 
Thus  the  only  essential  parameters  are  the  dimensionless  variances  ,  and 

C2<r>  anc*  higher  moments  characterizing  the  initial  data  of  the  bubble  size  and 

<jr2>  —  <  »  >2 

height  distributions.  Here  var(x)  =  - <x> -  Next  wc  cons^er  the  m°de'  >n 

an  infinite  strip,  jr0  2£  x  ss  x\.  The  ratio  of  mean  bubble  size  to  interval  width, 
xi  -  xo  is  a  dimensionless  measure  of  the  distance  of  the  data  from  the  final  equili¬ 
brium  state  of  a  single  bubble.  If  the  strip  is  finite,  zoS;S  z  i,  then  another  dimen¬ 
sionless  parameter  is  introduced  which  measures  the  time  at  which  the  interface  exits 
as  a  fraction  of  the  time  to  achieve  the  single  bubble  equilibrium  state. 

The  essential  phenomena  which  this  model  aims  to  capture  is  the  increase  in 
average  bubble  size,  and  the  consequent  increase  in  average  bubble  velocity,  resulting 
from  the  process  of  amalgamation.  In  order  to  study  the  amalgamation  quantita¬ 
tively,  we  introduce  a  growth  rate 

»  _  d  In  r 

k  ~  ~3T  • 

where  r  is  a  typical  bubble  radius.  From  (3),  we  see  that 


1 

v  St 


where  v  =  z  is  a  typical  bubble  velocity. 


The  mean  radius  of  the  bubbles  increases  in  time  because  the  model  is  based  on 
bubble  mergers  and  has  no  mechanism  by  which  a  bubble  can  split  into  smaller  bub¬ 
bles.  As  a  result,  the  possible  patterns  of  growth  of  the  ensemble  of  bubbles  can  be 
conveniently  discussed  with  reference  to  two  extreme  cases;  uniform  growth  and 
runaway  growth.  Uniform  growth  is  the  analog  of  a  steady  state  and  is  characterized 
by  a  translation  of  the  ensemble  distribution  in  In  r  with  time.  Runaway  growth  is 
the  approach  to  a  fixed  point  consisting  of  one  large  bubble.  Each  of  these  cases  may 
be  characterized  by  the  dynamics  of  the  support  of  the  ensemble  of  bubble  radii  in 
In  r  space.  In  the  case  of  uniform  growth,  the  support  of  the  bubble  distribution  in 
In  r  space  has  a  bounded  length,  which  is  0(1)  in  time.  At  the  other  extreme  is  the 
case  of  runaw-ay  growth,  in  which  the  support  of  the  bubble  distribution  shows 
marked  increase  in  length  as  a  function  of  time  because  the  radius  of  the  large  dom¬ 
inant  bubble  grows  at  a  greater  rate  than  the  mean  radius.  To  make  these  ideas 
quantitative,  we  define  the  length  of  the  support  to  be  In  rmax  -  In  r  where  rmax  is 
the  radius  of  the  largest  bubble  in  the  ensemble  and  r  is  again  a  typical  bubble 
radius. 

We  also  distinguish  two  time  periods  in  the  evolution  of  the  bubble  ensemble. 
We  suppose  initially  that  the  bubbles  are  generated  by  a  random  process,  and  thus 
that  the  distribution  of  a  given  bubble  (i.e.  its  radius  and  height,  considered  as  ran¬ 
dom  variables)  is  independent  of  its  neighbor.  In  other  words  adjacent  bubbles  are 
uncorrelated,  and  so  the  truncated  two  point  correlation  function  is  initially  zero. 
Under  this  hypothesis,  we  give  an  approximate  analysis  which  shows  that  the  regime 
of  uniform  growth  is  unstable  and  that  runaway  growth  occurs.  We  also  present 
numerical  evidence  which  shows  this  trend  very  clearly. 

However,  correlations  between  the  sizes  of  neighboring  bubbles  do  arise  dynam¬ 
ically  even  if  they  are  missing  initially.  In  fact  a  large  bubble  will  grow  rapidly  in 
radius  by  the  merger  process  when  placed  in  a  field  of  much  smaller  bubbles,  but 
rather  slowly  when  placed  in  a  field  of  only  slightly  smaller  bubbles.  Thus  after  an 
initial  period,  large  bubbles  are  more  likely  to  have  large  neighbors:  they  expand 
until  this  occurs  to  limit  or  retard  their  growth.  In  this  case  the  runaway  behavior  is 
self-limiting  and  a  regime  of  uniform  growth  appears  to  be  stable.  We  present 
numerical  evidence  for  the  stabilization  of  bubble  runaway  through  the  mechanism  of 
dynamically  generated  nearest  neighbor  correlations. 
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First  we  examine  an  initial  period  of  uncorrelated  nearest  neighbor  distributions. 
To  better  understand  the  behavior  of  the  two  regimes  of  uniform  and  runaway 

growth,  we  derive  an  approximate  formula  for  the  growth  rates  X  and  Xmax  =  max  . 

r  max 


Consider  a  local  region  of  space  in  which  a  single  large  bubble  is  growing  in  a 
uniform  background  field  of  small  bubbles  of  radius  r.  We  also  assume  a  small  dis¬ 
tribution  of  heights  for  the  background  bubbles.  Because  the  velocity  increases  with 
radius,  the  large  bubble  will  eventually  get  ahead  of  the  background  bubbles.  The 
height  of  the  large  bubble  is  then  fixed  at  a  height  C2r  above  the  height  of  the  small 
bubbles.  This  constraint,  combined  with  the  conservation  of  mass  in  the  merger  pro¬ 
cess,  forces  the  rapid  growth  of  the  large  bubble  to  be  primarily  in  the  lateral  direc¬ 
tion.  The  rate  at  which  bubble  mass  is  created  in  a  large  bubble  is  2rmaxciVgrraax  . 

Normalizing  for  the  effect  of  motion  of  the  small  bubbles  gives  a  corrected  mass  rate 
l  l 

of  2rmaxci[(grraax)7  -  (gr)7].  Thus  the  time  rate  of  change  2 rmax  of  the  diameter  is 
obtained  by  dividing  the  latter  expression  by  the  difference  in  height  between  the 
small  and  large  bubbles,  cir.  Thus 


Xmax  =  T225-  =  77''~1[(S'raax)T  ~  (gr)7]  . 

r  max  C  2 

When  rraax  »  r  ,  equation  (3)  implies  the  approximate  formula 

X  —  C ^  rmax  g 
*raax  r  -  v 

c2  r  vmax 

We  note  that  the  growth  rate  Xraax  goes  to  <*  with  the  aspect  ratio,  if  either  rmax  or  r 
is  held  fixed.  Hence  the  name  "runaway  mode"  to  discuss  the  case  in  which  this 
occurs.  It  is  convenient  to  consider  In  r  rather  than  r  as  the  basic  variable  in  terms  of 
which  the  radial  distribution  is  specified.  Thus  we  define  (Bln  r)max  through  the  for¬ 
mula 


t max  =  exp(ln  r  +  (Bln  r)max) 


and  then 


Araax  c 2  ^  r  ' 


exp( 


(Bln  r) 


max 


)  - 1 
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To  complete  our  semi-quantitative  analysis  of  the  stability  of  the  bubble  ensem¬ 
ble,  we  derive  a  formula  for  the  growth  rate  X.  Consider  two  adjacent  bubbles  hav¬ 
ing  radii  r+  and  r_  respectively.  For  simplicity,  we  suppose  that  initially  the  bubbles 
have  equal  heights.  We  determine  the  time  to  merge,  tm,  for  this  pair,  and  postulate 
that  this  time  is  the  doubling  time  of  the  mean  radius  r.  On  this  basis  we  determine 
the  growth  rate  X  in  the  following  equations 

2r  —  e  "  r 


or 


X  = 


In  2 


Moreover,  tm  is  determined  from  the  equations  of  motion  of  the  model,  (3)-(4),  to  be 

i-i 


We  choose  r *  as 


Then 


and  taking  ratios. 


or 


Cl 


1  1 
(gr+)T  -  (gr_)T 


r_  . 


r±  -  expdn  r  ±  r'> 


X-flnJ  (rt-rt)r:> 


(  rmax  'XT  _  , 

Lmax  _  1  r  -  Vr* _ 

*  *n  ^  r  ‘  x+  \y  /  r-  \y 


(^)T  *  (^) 


(  (Sin  r)max  \  , 

exp^ - j - )  -  1 


wmax 


-  1  rvr7  Sin  r\ 

_  In  2  exH  4  ) 


p(Sf>)  - . 


(5) 


From  these  formulas,  we  see  that  two  essential  parameters  which  describe  the 
input  statistics  have  a  critical  bearing  on  the  stability  (uniform  or  runaway  growth)  of 
the  bubble  merger  process.  These  parameters  are  the  variance  and  the  maximum 
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deviation,  Sin  r  and  (Sin  r)m „  of  the  In  r  distribution.  Similarly  the  variance  and 
maximum  deviation  of  the  height  /  radius  distribution  are  essential  parameters.  In 
addition  to  acting  as  above  as  a  cause  of  instability,  the  height  variation  can,  at  least 
initially,  contribute  to  stability.  The  mechanism  is  as  follows.  If  small  bubbles  are 
placed  ahead  of  and  adjacent  to  large  ones,  they  will  retard  temporarily  the  merger 
process  of  those  larger  bubbles. 

X, 


Analysis  of  (5)  shows  stability  (— ^ —  <  1)  or  uniform  growth  for 

(Sin  r)max 


ThT 


~  1  . 


(6) 


However  (6)  cannot  be  achieved,  except  for  the  identically  uniform  distribution  of 
bubbles  in  which  each  has  the  same  radius  and  height.  The  formulas  (5)-(6)  are 
intended  to  be  interpreted  qualitatively  rather  than  literally  and  thus  we  conclude  that 
control  of  rmax  relative  to  r  will  be  important  in  reducing  the  growth  rate  of  the  insta¬ 
bility,  and  in  delaying  the  onset  of  the  runaway  mode. 

If  the  initial  distribution  of  In  r  or  h  has  unbounded  support  (in  the  large  r  or  h 
directions)  then  outlying  large  bubbles  will  necessarily  occur  and  the  infinite  volume 
limit  is  necessarily  unstable  with  an  unbounded  growth  rate.  In  particular  it  would 
appear  that  the  runaway  period  would  not  be  self  limiting  a^d  that  the  subsequent 
period  of  uniform  growth  would  not  occur.  Our  computations  have  not  been  suffi¬ 
ciently  extepsive  to  determine  whether  uniform  growth  will  always  occur  for  the  case 
of  an  initial  distribution  with  bounded  support. 

Next  we  discuss  the  relation  between  the  infinite  volume  limit  and  a  considera¬ 
tion  of  the  model  in  a  bounded  strip  xoSx  s  jtj.  The  effect  of  the  finite  volume  is 
in  general  terms  to  convert  absolute  statements  into  probabilities,  so  that  for  example 
it  would  seem  that  for  a  finite  volume  and  an  initial  distribution  with  unbounded  sup¬ 
port,  a  finite  (non-zero)  probability  exists  that  the  self-limiting  uniform  growth 
region  will  not  occur. 

We  illustrate  these  ideas  by  the  analysis  of  a  series  of  computations.  Numeri¬ 
cally  we  observe  a  tendency  for  the  runaway  mode  to  be  self-limiting.  Large  bubbles 
expand  until  they  acquire  large  neighbors,  after  which  uniform  growth  occurs.  The 
details  of  this  phenomena  depend  sensitively  on  the  large  bubble  tail  of  the  distribu¬ 
tion,  and  so  we  do  not  yet  have  an  accurate  predictive  capability  for  this  region. 
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We  analyze  two  computational  solutions  of  the  bubble  model,  labeled  H  (homo¬ 
geneous)  and  L  (with  a  length  scale),  in  some  detail.  The  initial  distribution  of 
heights  h  and  radii  r  was  chosed  to  be  uniform  over  the  interval 

h  $  [ho  -  bh,  ho  +  bh],  r  €  [ro  -  br,  ro  +  Sr]  . 

Care  was  taken  in  generating  the  initial  random  interface  to  ensure  that  there  were  no 
sequential  correlations  among  the  bubbles.  The  units  for  both  runs  were  chosen  such 
that  the  initial  average  radius  was  <r>  =  1  ,  C2  =  1  ,  and  the  initial  mass  weighted 
average  velocity  was  <i>  =  1  .  The  initial  average  height  was  chosen  to  be 
<h>  —  0.  Both  runs  had  noisy  height  variation  bh/ro  —  1  .  Run  H  had  noisy  radial 
variation  8r/ro  =  1  ,  while  run  L  had  small  radial  variation  8r/ro  =  0.1  .  To  minim¬ 
ize  statistical  fluctuations  the  computations  were  started  with  about  10,000  bubbles. 

The  run  H  shows  an  initial  runaway  region  followed  by  uniform  growth,  as 
described  above.  Run  L  shows  two  new  features  in  addition  to  these  general  trends. 
Because  of  the  small  initial  variation  in  the  bubble  radius,  the  bubble  radius  distribu¬ 
tion  has  a  "resonance"  behavior,  or  a  multimodal  nature,  clustered  about  the  values 
nro,  for  the  bubbles  formed  from  the  merger  of  n  initial  bubbles.  The  second  prom¬ 
inent  feature  of  this  run  is  slow  start  and  the  extended  duration  of  the  initial  runaway 
period,  and  the  delayed  onset  of  the  self-limiting  uniform  growth  behavior. 

In  Fig.  4.1H  and  Fig.  4.1L,  we  show  the  location  in  x,  t  space  of  a  sampling  of 
new  bubbles  as  they  are  created  by  the  merger  of  two  or  more  neighboring  bubbles. 
Initially  the  mergers  are  randomly  distributed.  As  time  progresses  the  width  of  the 
bubbles  increases.  The  radial  distribution  of  bubbles  in  In  r  space  changes  with  time. 
A  sequence  of  plots  at  successive  time  steps  is  shown  in  Fig.  4.2  .  At  first  the  sup¬ 
port  of  the  distribution  grows  and  then  the  shape  of  the  distribution  becomes  smeared 
out. 

Many  aspects  of  the  merger  process  can  be  understood  from  time  histories  of 
appropriate  quantities.  Fig.  4.3  is  a  plot  of  the  log  of  the  number  of  bubbles  vs. 
time.  After  an  initial  period  the  number  of  bubbles  decreases  exponentially.  The 
initial  period  is  relatively  large  when  there  is  a  small  variation  in  the  radius  and 
hence  the  velocities  of  the  bubbles.  The  mass  weighted  average  velocity  vs.  time  is 
shown  in  Fig.  4.4  .  After  an  initial  period  of  slow  increase  the  velocity  accelerates 
almost  uniformly.  This  is  an  important  consequence  of  the  merger  process.  Mergers 


cause  the  mean  bubble  radius  to  grow  which  results  in  larger  velocities.  Similar  pred¬ 
ictions  of  uniform  acceleration  have  been  stated  previously.1,814  With  the  simple 
mode!  used  here  we  are  able  to  study  the  statistics  of  a  large  number  of  bubbles.  We 
find  the  magnitude  of  the  acceleration  depends  sensitively  on  the  bubble  distribution. 
In  Fig.  4.5  the  log  of  the  minimum,  mean,  and  maximum  bubble  radius  vs.  time  is 
plotted.  This  shows  the  growth  of  the  bounds  on  the  support  of  the  radial  distribu¬ 
tion  function.  Of  particular  interest  for  the  analysis  of  runaway  growth  is 
ln(rmM/rmean)  vs.  time  shown  in  Fig.  4.6.  We  see  that  the  support  grows  rapidly  indi¬ 
cating  a  runaway  growth  but  then  saturates  resulting  in  a  more  uniform  growth.  In 
Fig.  4.7  we  plot  the  minimum,  mean,  and  maximum  height  vs  time.  We  call  atten¬ 
tion  to  the  fact  that  the  maximum  height  is  the  quantity  which  would  be  relevant  to 
the  break  up  of  a  shell  for  example  in  the  laser  fusion  problem,  and  that  the  max¬ 
imum  behaves  qualitatively  differently  from  the  mean.  Furthermore,  this  behavior 
depends  on  the  initial  bubble  distribution.  The  spread  in  height  in  units  of  the  mean 
radius  is  very  much  larger  due  to  the  runaway  growth  when  the  initial  variation  in 
radius  of  the  bubbles  is  small  than  for  an  initial  noisy  variation  in  the  bubble  radius. 
The  saturation  in  the  growth  of  the  support  of  the  bubble  distribution  and  the  transi¬ 
tion  to  a  uniform  growth  is  due  to  correlations  which  evolve  dynamically.  In  Fig.  4.8 
we  plot  the  radius-height  cross  correlation,  or  in  other  words  the  correlation  between 
the  radius  and  the  height  of  a  single  bubble.  The  large  correlation  is  as  expected  and 
displays  the  fact  that  large  bubbles  move  faster  and  get  ahead  pf  small  bubbles.  In 
Fig.  4.9  the  nearest  neighbor  correlation  in  the  radius  and  height  is  shown, 

2(ri-  <r>)(ri+i-<r>) 

^-2(r,-<r>)2 

This  indicates  that  bubbles  grow  rapidly  until  the  adjacent  bubble  is  larger  than  the 
mean  in  radius  or  higher  than  the  mean  in  height.  Thus  the  runaway  growth  is  a 
local  phenomenon  that  depends  on  a  single  large  bubble  in  an  approximately  uniform 
background.  For  a  large  region  this  local  growth  pattern  may  occur  at  many  widely 
separated  places  which  eventually  interact.  Another  way  of  seeing  these  correlations 


which  limit  the  runaway  growth  is  to  consider  the  bubbles  adjacent  to  the  largest  bub¬ 
bles.  We  use  as  the  criterion  for  selecting  the  large  bubbles, 

In  r  >  0.2  In (r mean)  +  0.8  ln(rmax)  . 

In  Fig.  4.10  we  plot  the  average  bubble  radius  adjacent  to  .  :e  large  bubbles  in  units 
of  the  mean  radius  and  the  average  height  of  adjacent  bubbles  below  large  bubbles  in 
units  of  the  merge  height  for  bubbles  of  mean  radius.  Again  this  plot  indicates  that 
large  bubbles  run  into  adjacent  bubbles  with  larger  than  mean  radius  and  above  mean 
height.  These  runs  imply  that  the  dynamic  correlations  have  an  important  effect  on 
the  merger  process  and  that  a  simple  qualitative  model  as  outlined  in  Section  III 
would  require  a  non-trivial  closure  hypothesis  to  capture  this  phenomena. 

In  summary,  the  initial  period  of  runaway  behavior  is  characterized  by  a  small 
but  exponentially  increasing  bubble  acceleration  and  a  relative  absence  of  nearest 
neighbor  correlations.  The  transition  to  uniform  growth  is  caused  by  the  dynamic 
development  of  nearest  neighbor  correlations  and  is  characterized  by  a  relatively  con¬ 
stant  but  large  acceleration. 

V.  NUMERICAL  DETERMINATION  OF  THE  MODEL  PARAMETERS 

A.  The  One  Body  Problem 

We  identify  three  regimes  for  tjie  time  development  of  the  rising  bubble  and  for 
the  falling  spike,  in  the  case  of  motion  in  an  infinite  strip,  so  that  boundary  effects 
can  be  neglected.  The  initial  time  for  the  spike  and  bubble  motion  is  set  by  the  initial 
amplitude  of  a  sinusoidal  surface.  The  perturbation  of  the  initially  flat  interface  is 
small,  and  the  time  evolution  of  the  surface  can  be  described  by  the  linearized  Euler 
equations.  (See  the  Appendix.) 

The  isothermal  equilibrium  situation  consists  of  a  flat  interface  with  exponen¬ 
tially  stratified  atmospheres  above  and  below  the  interface.  We  consider  perturba¬ 
tions  of  the  unstable  equilibrium  in  which  the  heavy  fluid  lies  below  the  light  fluid 
and  gravity  is  directed  upwards. 

The  compressible  Rayleigh-Taylor  problem  depends  on  three  dimensionless 

Pi 

parameters:  the  density  ratio  D  -  — — ,  where  p*  is  the  density  of  the  heavy  gas  just 


below  the  interface  and  pa  is  the  density  of  the  light  gas  just  above  the  interface,  the 
polytropic  gas  constant  y  (here  we  set  ya  =  7*)  or  other  information  to  set  the  equa¬ 
tion  of  state  for  the  heavy  and  light  fluids,  and  a  parameter  M  defining  the  ratio  of  a 
gravitational  time  scale  to  a  sound  speed  time  scale.  ( M  defines  a  dimensionless 

compressibility19.)  We  take  M 2  =  -^j-,  where  X  is  the  wavelength  of  the  interface 

perturbation  and  q,  is  the  sound  speed  in  the  unperturbed  heavy  fluid  In  place  of  the 

density  ratio  D,  the  Atwood  ratio  A  =  —  mav  be  used  as  a  dimensionless 

Pi  +  Pa 

parameter.  Observe  that  the  mass  of  the  heavy  fluid  per  unit  cross  section  is  finite, 
even  though  the  strip  it  is  located  in  is  semi-infinite,  due  to  the  exponential  stratifica¬ 
tion  of  the  fluid  densities.  In  fact  this  mass  density  is 

at 

fi  =  Pi  ^exp(-Py)  dy  , 

where  B  =  . 

Cb 2 

The  initial  growth  at  small  time  of  both  the  spike  and  bubble  is  exponential. 
The  formulation  of  the  linear  stability  theory,  which  governs  the  small  time  behavior, 
was  given  by  Bernstein  and  Book1*  for  the  infinite  strip.  These  results  are  given  in 
the  Appendix,  and  are  extended  to  apply  to  a  finite  strip  with  reflecting  boundary 
conditions  at  top  and  bottom.  They  are  also  extended  to  apply  to  a  general  equation 

*  CfrG 

of  state.  The  dependence  of  the  dimensionless  exponential  growth  rate  constant  — — 

O 

on  the  dimensionless  parameters  A  and  M  is  shown  in  Fig.  5.1  for  y  ~  1.4.  In  Fig. 
5.1,  the  dimensionless  growth  rate  contours  are  equally  spaced,  with  the  difference 
between  two  adjacent  contours  being  .793.  A  larger  value  of  M2  or  A  gives  rise  to  a 
larger  dimensionless  growth  rate,  since  gravity  increases  with  M2  and  renormalized 
gravity  increases  with  A.  At  intermediate  times,  the  spike  enters  a  period  of  free 
fall.  The  bubble  may  have  an  intermediate  period  of  free  rise,  and  if  it  does  so  it  is 
or  may  be  Rayleigh-Taylor  unstable  and  is  therefore  potentially  unstable  to  splitting 
through  the  development  of  a  spike  at  its  center.  At  late  times  both  the  bubble  and 
the  spike  approach  an  asymptotic  constant  velocity  due  to  the  form  drag  and  the  finite 
total  mass  supported  above  the  bubble  and  spike.  Although  not  present  in  the  model 
equations  discussed  here,  an  experimentally  more  realistic  late  time  behavior  for  the 


spike  would  be  a  breakup  into  droplets  which  then  also  achieve  a  terminal  velocity 
due  to  form  drag.  The  regimes  of  free  fall  and  of  the  approach  to  an  asymptotic 
velocity  for  the  spike  can  be  modeled  by  a  particle  falling  under  the  influence  of 
gravity  g  with  a  drag  coefficient  b  and  limiting  velocity  vterm.  The  model  problem 

v  =  -  b(v  -  vterm)  (7) 

has  the  first  integral 

v  =  vterm  +  (v0  -  vterm)  exp (~b(t  -  t0 ))  (8) 

where  b  is  the  drag  coefficient,  vo  is  the  velocity  at  time  to  and  vterm  is  the  terminal 
spike  velocity.  Actually,  equation  (7)  is  the  expansion  of  the  dynamic  equation  in  a 
neighborhood  of  the  terminal  velocity  keeping  only  the  leading  term  (the  term  linear 
in  the  velocity).  Therefore  equation  (8)  should  be  valid  only  in  the  terminal  velocity 
region.  To  understand  the  behavior  of  the  velocity  over  a  larger  domain  (including 
the  linear  and  free  fall  regions),  the  contribution  from  the  higher  order  terms  should 
also  be  included.  To  the  third  order  in  velocity,  the  dynamic  equation  can  be  written 
as 

V  =  fllV  +  d2V2  +  tf3V3.  (9) 

One  can  see  by  inspection  that  equation  (9)  has  qualitatively  the  desired  form.  Obvi¬ 
ously,  it  has  a  pritical  point  v  =  0.  By  suitable  choices  of  a i,  a i  and  <23,  it  will  have 
another  critical, point  |v|  =  |vterra|  >  0.  In  general,  a  third  order  algebraic  equation 
should  have  three  solutions.  But,  as  we  will  see,  one  of  the  solutions  is  not  physical. 
The  critical  velocity  v  =  0  corresponds  to  the  velocity  for  a  perturbation  with  infini¬ 
tesimal  initial  amplitude  at  the  interface.  The  system  will  leave  this  critical  point 
exponentially  and  it  approaches  the  v  =  critical  point  exponentially  also.  We 
now  develop  these  facts  systematically,  and  also  determine  how  the  terminal  velocity, 
renormalized  gravity  in  the  free  fall  region  and  drag  coefficient  are  related  to  the 
coefficients  of  equation  (9). 

In  the  earlier  stage  of  the  development,  one  may  neglect  the  terms  proportional 
to  v2  and  v3  in  equation  (9)  due  to  the  small  amplitude  of  the  velocity.  Then  equa¬ 
tion  (9)  can  be  approximated  as 
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This  equation  has  the  solution 

v  =  v0  exp(ai(r  -  t0  )).  (11) 

We  conclude  a\  =  a,  by  comparison  to  the  results  from  the  linear  theory.  (See  the 

a2 

Appendix.)  By  a  dimensional  argument,  ai  must  have  the  form  <22  =  —  2k —  and 

5 

CT3 

a 3  =  -  £— > r,  where  tc  and  £  are  dimensionless  constants  and  g  is  (the  unrenormal- 
g 

ized)  gravity. 

The  concavity  of  the  curve  representing  velocity  versus  time  is  determined  by 
the  derivative  of  equation  (9),  i.e., 

4"r  =  (CT  -  4k-^-v  -  3£-^j-v2)v.  (12) 

at  g  g 


Then  the  curve  at  velocity  v  is  concave  up  or  down,  depending  on  whether  the 
derivative  of  the  acceleration  is  positive  or  negative.  At  the  inflection  point  the 
derivative  of  the  acceleration  vanishes.  If  we  define  a  velocity  Vf,  determined  by  the 
equation 

a  _  4K£-Vf  -  3^-p-Vf2  =  0,  (13) 

then  the  linear  theory  region  corresponds  to  the  range  where  v  «  Vf,  the  terminal 
velocity  region  corresponds  to  the  range  where  v  »  Vf,  and  the  free  fall  region 
corresponds  to  the  range  where  v  ~  Vf. 

The  terminal  velocity  vterra  is  determined  by  setting  the  acceleration  in  equation 
(9)  to  zero,  i.e. 

(Tv  -  2k— v2  -  =  0.  (14) 

8  g2 

Obviously,  v  =  0  is  not  a  solution  for  the  terminal  velocity.  The  other  two  possible 
solutions  are 


1 

V-  =  -^(K  +  (K2  +  0T] 


(15) 


1 

V+  =  -zf[K  '  (*2  +  *)T}' 


(16) 


and 


Physically,  the  system  starts  with  a  perturbation  at  the  interface.  The  initial  velocities 
of  the  spike  and  bubble  are  determined  by  the  linear  theory.  The  terminal  velocity 
must  have  the  same  sign  as  the  initial  velocity.  If  both  solutions  vs  satisfy  this  sign 
condition,  then  the  solution  with  the  smaller  magnitude  is  the  terminal  velocity,  since 
the  system  starts  at  the  critical  point  v  ■  0  and  will  meet  the  smaller  amplitude  solu¬ 
tion  first.  Let  us  consider  the  spike  first.  We  choose  0  <  g.  For  0  <  £,  only  v+  has 
the  same  sign  as  the  initial  spike  velocity  (0  <  v+).  For  £  <  0  and  k  <  0,  both  v~ 
and  v“  are  negative  and  neither  satisfies  the  sign  condition.  For  £  <  0  and  0  <  k, 
both  v+  and  v“  satisfy  the  sign  requirement,  but  v*  has  the  smaller  amplitude 
Therefore  v+,  if  it  is  positive,  is  the  only  physical  solution  for  the  terminal  velocity 
of  the  spike,  while  v~  is  always  an  unphysical  solution.  Similarly,  v“,  if  it  is  nega¬ 
tive,  is  the  only  physical  solution  for  the  terminal  velocity  of  the  bubble  and  v*  is  an 
unphysical  solution.  Furthermore,  for  sufficiently  negative  values  of  <;  (£  <  -k2), 
both  v*  and  v~  become  complex  quantities.  In  this  case  the  terminal  region  does  not 
exist  for  the  model  system  (9)  in  that  parameter  range.  For  the  same  reason  equa¬ 
tion  (13)  also  has  at  most  one  physical  solution.  For  the  spike 

_  l 

vf  =  3^-[2k  -  (4k2  +  35)T],  (17) 

and  for  the  bubble  \ 

_  1 

vf  =  3^[2k  +  (4k2 -f  3€)7].  (18) 

Moreover  renormalized  gravity  in  the  free  fall  region  is  given  by 

gR  =  I  <rvf  -  2k~v f2  -  t-fjV  |.  (19) 

®  a 

We  rewrite  equation  (9)  in  terms  of  the  terminal  velocity.  Then  the  coefficient 
of  the  linear  term  is  the  drag  coefficient,  i.e. 

UP"  ~  b(v  ~  vterm)  +  0[(v  —  Vtenn)2] 


b  —  (<T  4k  v^rm  3£— yVt2erm) 

®  8 


where 


Here  the  plus  sign  is  for  the  bubble  and  the  minus  sign  is  for  the  spike. 


The  solution  for  equation  (9)  is 

t  —  tQ  =  ■  £5*  [  _  ln(-j-4-)  +  —  - - - 

£<X3  V  V  |vo|  V  (v  —  V’ 

I  _  I 

1 


V,  -  V 


r’n(-r - : 

)  |v'o  -  V 


+ 


Trln(- 


V,  -  V 


V  (v  -  O  | v 0  -  V  j 


T=T>]- 


(21) 


where  v,  =  v(r)  is  the  velocity  at  time  t,  and  vT  and  v“  are  given  by  equations  (15) 
and  (16).  This  is  the  formula  we  use  to  fit  the  velocities  of  the  bubble  and  spike 
simultaneously  over  all  three  different  (linear,  free  fall  and  terminal)  regions.  There 
are  two  free  parameters  in  (21)  which  can  be  adjusted  to  fit  the  computed  curve  of 
spike  velocity  vs.  time,  as  shown  in  Fig.  5.2.  The  computed  spike  and  bubble  veloci¬ 
ties  were  smoothed  by  averaging  before  plotting  to  eliminate  oscillations  of  a  numeri¬ 
cal  nature  in  the  data. 

In  Fig.  5.3,  we  illustrate  the  spike  dynamics  in  three  (linear,  free  fall  and  termi¬ 
nal)  regions.  Renormalized  gravity  corresponds  to  the  maximum  acceleration  over 
the  whole  time  domain.  This  is  a  consequence  of  our  discussion  concerning  the  con¬ 
cavity  of  the  velocity  curve.  When  the  system  goes  from  the  linear  region  to  free 
fall,  the  deviation  of  the  linear  theory  from  the  results  for  the  full  Euler  equation 
becomes  substantial.  Likewise  before  the  system  enters  the  terminal  region,  the 
results  from  the  asymptotic  behavior  (Eq.  (7))  give  a  more  rapid  deceleration  than 
the  true  value. 

In  Fig.  5.4,  we  present  the  dependence  of  the  gravity  renormalization  factor 
gt(/g  defined  by  equation  ( i9:  on  the  density  ratio  D.  The  data  (from  reference  20) 
for  two  dimension  incompressible  ( M 2  =  0),  one  component  ( A  =  1)  flow  are  also 
shown.  The  incompressible  data  may  be  affected  by  the  large  initial  amplitude  (0.5) 
used  in  the  cited  computation.  Our  computational  data  show  that  the  renormalization 
factor  depends  linearly  on  the  Atwood  ratio  A  and  g*  ~  .5 Ag  for  the  bubble,  while 
gK  ~  Ag  for  the  spike. 

At  equilibrium,  the  right  side  of  equation  (15)  vanishes,  and  in  this  case  we 
have  an  expression  for  the  drag  d  =  Z.i+,ftg  where  Lt,+S  is  the  combined  width  of  the 
bubble  and  spike,  ft  is  the  mass  density  given  at  the  beginning  of  this  section.  It  is 
convenient  to  introduce  the  dimensionless  drag  coefficient  of  the  spike  following 
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reference 


20 


CD  =  2 


Pav?errt/-j 


(22) 


where  Lj  is  the  width  of  the  spike  tip.  Substitution  of  the  expression  for  d  into  above 
equation  yields 


C„  =  2  % 

Ls  pvferm 


(23) 


The  computed  drag  coefficient  Co  as  a  function  of  A  and  M  is  plotted  in  Fig.  5.5  for 
•y  =  1.4. 

The  terminal  Mach  number,  vterm/cb,  is  plotted  as  a  function  of  the  Atwood  ratio 
in  Fig.  5.6.  Here  c*,  is  the  sound  speed  of  the  heavy  material  just  below  the  tip  of  the 
spike  (Fig.  5.6a)  or  the  bottom  of  the  bubble  (Fig.  5.6b).  We  note  that  the  terminal 
Mach  number  has  a  significant  dependence  on  the  compressibility  M2.  The  data  for 
D  =  10  and  100  are  not  included  in  the  plot  because  these  systems  did  not  reach  the 
terminal  velocity  region  at  the  final  computational  step.  At  the  termination  of  these 
runs,  the  dimensionless  heights  which  are  the  ratio  of  the  amplitude  of  bubble  or 
spike  to  the  wave  length  are  in  the  range  from  1  to  2.5,  but  the  systems  are  still  fully 

in  the  free  fall  regime.  The  spike  and  the  bubble  velocities  do  become  supersonic 

\ 

(relative  to  the  sound  speed  of  the  heavy  gas)  before  the  final  computation  step  in  the 
D  =  10  and  100,  M2  =  2  cases.  For  D  =  100,  M2  =  2  the  spike  Mach  number 
reaches  2.7,  the  bubble  Mach  number  reaches  2.3  and  they  are  still  growing.  How¬ 
ever  in  all  cases  considered  here,  the  bubble  and  spike  velocities  are  subsonic  relative 
to  the  light  gas.  If  the  light  fluid  is  a  vacuum  (D  =  «,  A  =  1)  then  the  spike  is 
expected  to  remain  in  free  fall.  This  behavior  is  consistent  with  the  behavior  of  an 
incompressible  fluid.21 

A  dimensional  analysis  of  the  Euler  equation  shows  that  the  terminal  velocity 
vterm  of  the  bubble  scales  according  to  equation  (3) 


A 

Krai  =  Q  (gr)T, 


(24) 


where  c\  is  a  function  of  the  fundamental  dimensionless  parameters  D,  M  and  y 
alone  and  r  is  half  wave  length,  i.e.  the  radius  or  half  the  width  of  the  single  mode 
computation.  Fig.  5.7  illustrates  the  dependence  of  the  constant  c\  on  D  and  A/2.  As 
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noted  above,  the  systems  with  high  density  ratios  have  not  reached  the  terminal  velo¬ 
city  region  at  the  end  of  the  computation.  Therefore  the  terminal  velocity  can  not  be 
determined  in  these  cases.  However,  the  velocity  at  the  last  step  of  the  computation 
provides  us  with  a  lower  bound  for  cj  in  these  systems.  For  D  =  5,  A/2  =  .5, 
Ci  >  .36,  for  D  -  10,  Af 2  =  .5,  ci  >  .46  and  for  D  =  5,  A/2  =  2,  ci  >  .47.  The 
data  for  a  two  dimensional  (from  reference  21  and  22)  incompressible  (A/2  =  0),  one 
component  ( A  -  1)  fluid  is  also  shown  in  Fig.  5.7.  There  is  an  expectation  that 

l 

ci  ~  const  4^  In  fact  for  both  the  linear  and  free  fall  regimes,  gravity  occurs  in  the 

form  Ag.  If  this  remains  the  case  for  the  terminal  velocity  regime  in  equation  (24), 

l 

then  ci  would  be  proportional  to  A^ .  Our  data  indicates  that  ci  has  a  significant 
dependence  on  A/2.  Even  allowing  for  an  A/2  dependence  in  the  constant  of  propor¬ 
tionality,  the  data  does  not  fit  this  relation  particularly  well,  especially  with  the  inclu¬ 
sion  of  the  lower  bounds  for  systems  not  in  the  terminal  velocity  region. 

The  interaction  of  the  spike  and  bubble  with  waves  reflected  from  the  boun¬ 
daries  leads  to  a  slowing  down  of  their  motion,  and,  for  motion  in  a  bounded  region, 
defines  a  further  time  regime  for  the  fluid  motion. 

We  have  made  a  series  of  computations  sampling  the  development  of  the  single 
finger  as  a  function  of  D.  Figs.  5.8-5.10  present  plots  of  the  interface,  density  con¬ 
tours,  and  pressure  contours  for  D  =  1.5,  10,  and  100  and  A/2  =  .5.  Fig.  5.11-12 
presents  plots  of  the  interface  for  D  =  2  and  100  and  Af2  =  2.  Note  that  the  flow  is 
subsonic  in  Fig.  5.8  and  supersonic  relative  to  the  heavy  gas  in  the  final  frame  of  Fig. 
5.12 

Because  there  are  two  dimensionless  parameters  in  the  formulation  of  the  prob¬ 


lem,  D  and  Af2  (with  the  equation  of  state  held  fixed),  there  should  also  be  two  quali¬ 
tative  features  in  morphology  of  the  fingers  in  the  deeply  nonlinear  regime,  which 
change  in  their  relative  importance  as  D  and  Af2  vary.  The  effect  of  increasing  D  is 
seen  in  a  trend  toward  a  thinner  spike  and  less  roll  up  shed  off  the  edge  of  its  tip. 
The  effect  of  increasing  A/2  leads  to  the  deposit  of  a  strip  or  trailing  filament  of 
material  shed  off  the  edge  of  the  tip  in  contrast  to  the  highly  wound  up  vortices 
which  are  observed  in  the  incompressible  case.  The  increase  of  A/2  leads  to  a 


a 


dramatic  increase  in  terminal  bubble  and  spike  Mach  numbers.  For  large  values  of 
Af2,  a  bow  shock  forms  in  the  front  of  the  bubble  and  a  complex  system  of  shocks 


mmm 


form  in  the  stem  of  the  spike.  The  present  calculations  are  only  moderately  super¬ 
sonic  and  so  the  shocks  are  weak.  However,  due  to  the  large  value  of  gravity  and  the 
exponential  stratification  of  the  pressure,  the  pressure  changes  very  rapidly  in  a  few 
mesh  layers  even  in  the  absence  of  shock  waves,  so  that  these  weak  shocks  are  not 
clearly  resolved  in  the  present  calculations. 

B.  The  Two  Body  Problem 

In  Fig.  5.13-14  we  show  a  series  of  frames  which  give  the  time  development  of 
a  bubble  merger.  The  parameters  for  this  solution  are  D  =  2,  10  and  M2  =  2  .  A 
detailed  study  of  the  two  body  problem  would  be  an  interesting  topic  for  future 
work,  and  is  required  to  complete  the  determination  of  the  parameters  in  the  Sharp- 
Wheeler  model  and  their  dependence  on  the  dimensionless  parameters  D  and  M2  of 
the  compressible  Rayleigh-Taylor  problem  from  first  principles.  In  any  case  we 
believe  that  the  preliminary  result  reported  here  indicates  that  this  goal  is  feasible, 
using  the  front  tracking  method. 

C.  Heterogeneity 

Heterogeneity  can  have  an  important  effect  on  solutions.  This  is  evident  from 
6ur  study7  of  the  incompressible  Rayleigh-Taylor  problem.  A  detailed  analysis  of 
this  topic  is  postponed  to  a  subsequent  paper,  but  for  the  present  we  observe  that 
heterogeneity  in  the  form  of  turbulence  in  the  flow  field  could  lead  to  pinch  off  of 
the  spike  or  of  the  trailing  edge  of  the  spike  roll  up,  as  well  as  the  break  up  of  the 
spike  into  droplets.  Any  of  these  phenomena  would  effect  the  drag  coefficient  and 
thus  the  terminal  velocity  of  the  spike.  We  expect  that  renormalized  gravity,  which 
characterizes  the  free  fall  regime,  should  be  less  sensitive  to  heterogeneity. 

These  preliminary  results  also  suggest  additional  phenomena,  such  as  bubble 
pinch  off,  true  multiphase  flow  and  a  possible  variety  of  multiphase  flow  regimes. 
The  bubble  mergers  presented  here  actually  suggest  the  regime  of  slug  flow,  which  is 
recognized  as  a  regime  in  the  multiphase  flow  of  fluids  in  pipes.  In  addition  there  is 
good  reason  to  suppose  that  neighboring  bubbles  influence  one  another.  For  exam¬ 
ple,  two  large  and  neighboring  bubbles  might  have  a  collective  terminal  velocity 
which  is  larger  than  either  would  have  individually.  The  extent  and  importance  of 
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such  additional  phenomena  is  left  to  further  studies  as  is  the  task  of  modifying  as 
necessary  the  statistical  model  for  bubble  merger. 

VI.  VALIDATION  OF  SINGLE  FINGER  COMPUTATIONS 

In  the  small  amplitude  regime,  the  computed  solution  can  be  compared  to  the 
analytic  solution,  derived  in  the  Appendix,  of  the  linearized  equations,  where  the 
finger  is  considered  as  a  small  amplitude  perturbation  of  the  flat  interface.  At  any 
amplitude,  mesh  refinement  produces  a  validation  test  of  the  solution.  In  addition 
we  have  determined  that  the  solution  is  independent  of  various  numerical  parameters 
which  adjust  the  algorithm,  and  some  of  the  these  results  are  presented  here.  In  par¬ 
ticular,  we  study  the  effect  of  the  initial  amplitude  of  the  perturbation,  the  effect  of 
the  boundary  condition  at  the  top  and  bottom  of  the  computational  domain  and  the 
effect  of  the  remeshing  (redistribution)  of  points  along  the  interface. 

We  study  first  the  effect  of  mesh  refinement.  Obviously  the  finer  the  mesh 
which  is  used,  the  greater  is  the  level  of  detail  which  one  can  see.  Convergence 
occurs  in  the  sense  that  structures  which  occur  on  the  coarse  grid  are  duplicated  on 
the  fine  grid.  In  Fig.  6.1  we  illustrate  mesh  refinement  by  a  factor  of  4  in  each  linear 
dimension.  We  note  that  the  degree  of  resolution  we  obtain  could  not  be  achieved  on 
comparable  grids  without  the  use  of  interface  methods.  In  this  figure  D  =  10  and 
M2= 0.5.  As  we  see  in  Fig.  6.2,  the  positions  of  the  spike  and  the  bubble  are  hot 
very  sensitive  to  mesh  refinement.  The  reason  is  that  convergence  for  these  features 
of  the  solution  has  largely  occurred  on  the  coarse  grid,  and  so  the  mesh  refinement 
only  affects  the  detailed  structure  of  the  interface  and  in  particular  the  secondary 
instabilities  on  the  side  of  the  spike.  Secondary  structures,  which  occur  on  the  finest 
grid  only,  as  in  Fig.  5.9  and  Fig.  6.1  are  not  confirmed  by  the  validation  studies 
presented  here.  They  are  somewhat  sensitive  to  numerical  parameters  such  as  fre¬ 
quency  of  remeshing;  frequent  remeshing  of  the  interface  is  equivalent  to  smoothing 
and  tends  to  suppress  these  secondary  structures.  We  anticipate  that  the  secondary 
structures  also  depend  sensitively  on  surface  physics  effects,  such  as  surface  tension 
or  a  viscous  diffusion  layer. 

In  the  small  amplitude  region,  the  interface  is  smooth  and  nearly  flat.  In  that 
region,  the  spike  and  bubble  grow  exponentially.  We  define  a  dimensionless 
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amplitude,  a,  the  ratio  of  the  initial  perturbation  amplitude  at  the  interface  to  the 
wave  length.  In  order  to  test  carefully  the  computation  in  the  linear  regime,  we  start 
at  a  very  small  value  of  a,  a  =  0.001,  in  Fig.  6.3  and  show  the  comparison  of  the 
bubble  and  spike  velocity  and  the  position  for  the  linear  theory  and  the  full  Euler 
equations.  In  Fig.  6.4,  we  have  plotted  the  velocity  and  the  position  of  the  bubble 
starting  with  two  different  values  of  a,  a  =  .004  and  a  =  .01,  superimposed.  As  we 
can  see  from  the  figure,  the  end  of  the  linear  regime  and  the  behavior  of  the  system 
after  exit  from  the  linear  regime  are  independent  of  the  initial  starting  amplitude  pro¬ 
vided  that  this  amplitude  lies  in  the  linear  region.  In  the  study  of  the  nonlinear 
regime,  it  is  desirable  to  start  the  initial  amplitude  as  large  as  possible  to  minimize 
the  effects  of  the  reflected  wave  from  the  boundary,  which  we  will  discuss  in  Fig. 
6.6. 

Because  of  the  agreement  between  the  computed  solution  and  the  linear  theory 
and  between  the  computed  solutions  with  different  starting  amplitudes,  we  believe 
that  the  initial  amplitude  and  the  linear  ansatz  in  the  initialization  of  the  computed 
solution  have  no  effect  on  the  computation.  Furthermore,  we  see  that  the  linear 
theory  gives  good  agreement  up  to  amplitudes  a  -  0.015  and  so  these  values  of  ini¬ 
tial  a  are  used  in  the  further  studies  of  the  nonlinear  regime. 

In  the  small  amplitude  regime,  the  interface  is  smooth  and  nearly  flat.  There¬ 
fore  the  remeshing  of  the  interface  is  not  required,  and  its' only  effect  is  a  regulariza¬ 
tion  through  interpolation.  Thus  we  see  in  Fig.  6.5  that  remeshing  at  a  low  fre¬ 
quency  (once  per  15  or  once  per  100  time  steps)  has  no  influence  on  growth  rates, 
while  a  remeshing  every  time  step  has  the  expected  effect  of  lowering  the  growth  rate 
of  the  instability.  Furthermore  a  frequent  remeshing  suppresses  the  secondary  insta¬ 
bilities  on  the  surface  of  the  spike. 

For  small  times,  the  upper  and  lower  boundaries  have  no  effect  on  the  motion 
of  the  interface  due  to  the  finite  propagation  speed.  In  the  free  fall  regime,  we  have 
observed  that  the  effect  of  the  reflected  boundary  signals  is  not  very  strong.  In  the 
terminal  velocity  region,  reflected  boundary  signals  can  have  a  very  noticeable  effect 
on  the  interface.  In  Fig.  6.6,  we  compare  the  bubble  position  and  velocity  for  a 
domain  of  shape  1  x  4  and  1  x  8.  At  late  times  the  boundary  signals  do  effect  the 
terminal  velocity.  The  times  marked  by  an  arrow  show  the  first  arrival  of  a  signal 


from  the  nearest  boundary  and  the  first  arrival  of  a  reflected  signal  starting  initially 
at  the  interface.  The  first  of  these  times  marks  the  beginning  of  a  quantitative  diver¬ 
gence  between  these  two  curves  and  the  second  marks  the  beginning  of  a  qualitative 
difference,  with  the  boundary  reflected  signal  causing  the  velocity  to  begin  decreasing 
in  magnitude  in  the  shorter  region.  The  other  anomalies  in  the  velocity  vs  time  plots 
are  coarse  mesh  phenomena,  and  clearly  have  nothing  to  do  with  boundary  effects,  as 
can  be  seen  by  the  fact  that  they  occur  at  the  same  locations  for  the  two  runs  with  dif¬ 
ferent  boundary  locations.  To  further  emphasize  this  point  we  also  plot  in  this  figure 
a  finer  grid  computation  of  the  same  quantities  for  the  shorter  domain. 

At  large  time,  it  is  necessary  to  remesh  the  interface  due  to  the  large  scale  dis¬ 
tortions  in  the  geometry.  It  seems  that  remeshing  at  on  the  order  of  at  least  every  20 
time  steps  is  required.  If  the  interface  is  remeshed  more  frequently,  then  the  secon¬ 
dary  structure  of  the  surface  waves  on  the  side  of  the  spike  are  suppressed  but  the 
position  of  the  bubble  and  spike  is  not  greatly  affected.  For  the  extreme  case  of 
remeshing  every  time  step,  however,  the  growth  rates  of  the  bubble  and  spike  are 
affected,  as  we  have  seen. 
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APPENDIX:  THE  LINEARIZED  COMPRESSIBLE  RA  YLEIGH-TAYLOR 
THEORY 

We  derive  the  equations  governing  the  initialization  of  the  single-mode 
Rayleigh-Taylor  computations.  Although  our  derivation  is  given  explicitly  for  the 
single-mode  problem,  the  solutions  for  the  multiple-mode  Rayleigh-Taylor  computa¬ 
tions  can  be  obtained  by  the  superposition  of  the  single-mode  solutions,  since  the 
equations  are  linear  and  different  modes  do  not  couple  with  each  other  at  this  order. 
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Our  analysis  extends  results  of  ref.  19. 

The  two-dimensional  Euler  equations  for  a  compressible,  in  viscid,  gas  are 

4®-  +  =  o  (Al) 

dt  dx  dz 

teL  +  a(PM*  +  +  tee.  =  0  (A2) 

dt  dx  dz 

Hr  +  *zr+  *(pv^— f)  =  Pi  (A3) 

^■(P(yi:  +  *)]  +  •§; [pw(y93  +  i))  +  ^-[p»'(y^2  +  0]  -  p vi,  (A4) 

where  u  is  the  x  component  of  the  velocity  ,  v  is  the  z  component  of  the  velocity, 
q2  =  u2  4-  v2,  e  is  the  specific  internal  energy  and 

i  =  *  +  j  (A5) 

is  the  specific  enthalpy.  The  thermodynamic  quantities  pressure  and  density  are 

related  by  the  equation  of  state, 

e  -  e(p, p).  (A6) 

Suppose  we  have  a  perfect  flat  interface  between  the  light  gas  (above  a)  and  the 
heavy  gas  (below  b\.  This  isothermal  equilibrium  state  is  determined  from  the  fol¬ 
lowing  equations: 
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IT  =  pog 
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O 

(A8) 

eo  -  e(Po,Po). 

(A9) 

Here  the  subscript  zero  represents  the  solutions  of  the  isothermal  equilibrium  state, 
and  gravitation  is  taken  to  point  in  the  positive  z  direction. 

To  analyze  the  Rayleigh-Taylor  instability,  we  imagine  that  the  interface 
between  the  light  and  the  heavy  gas  is  perturbed  as  a  sine  wave  with  one  mode.  Since 
the  perturbation  is  small,  the  state  of  the  system  should  be  close  to  the  solution  of  the 
isothermal  equilibrium  state  on  the  short  time  scale.  Therefore  to  linearize  the  equa¬ 
tions  (A  1)  — (A4)  and  (A6),  we  write  each  function  as  its  zeroth  order  (equilibrium) 
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solution  and  a  first  order  correction. 


p  ~  po  +  8p 


(A10) 


P  ~  P0  +  hP 


(All) 


u  ~  hu 


(A12) 


v  ~  8v 


(A  13) 


e  ~  eQ  +  8e. 


(A  14) 


Here  we  have  used  the  results  that  the  zeroth  order  solutions  for  velocities  in  the  x 


and  z  direction,  uq  and  vq,  are  zero.  The  linearized  equations  are 


88v  ,  dp0 , 


,  _  VKJU  ,  _  u  u  V  ,  £  —  A 

+  po_ar +  po~  +  ~8v  - 0 


(A15) 


„  58m  ,  dhp  _  n 

po-rr  +  tt  -  0 


(A  16) 


„  asv  ,  dhp  _  > 
"0—  +  TT  -  s8p 


(A  17) 


d8p  ,  „  dhe  ,  „  .  38u  ,  .  dv  ,  .  dpp  c 

e°~dT  +  Po~TT  +  p0'°^7"  +  p0,°  d?  +  l°~dT8v 


dl’o  _ 

+  P0-37&''  -  Po^ov 


(A  18) 


dfop  _  ^££^8g_  de0  dhP 

dt  dpo  dt  6Pq  dt 


(A  19) 


By  using  equations  (A5),  (A8),  (A15)  and  (A19),  equation  (A18)  can  be  rewritten  as 


d*Q  Po.  38o  ,  .  3*0  38P  _  .  c 
(poap^  ~  p0^~  _  poghv- 


dP  0  dt 


(A20) 


Since  all  the  coefficients  are  the  functions  of  z,  we  can  write  each  function  as 


8/,  =  exp(ar  +  ikx)hfi(z) , 


(A21) 


for  8/,  =  8p,8P,8u,8v.  Here  hf]  is  a  function  only  of  z.  Each  value  of  a  and  k 
corresponds  a  single  mode.  By  using  equation  (A21),  the  linearized  equations 


become 


o8p  +  ikpohu  +  po~^~  +  4jp-8v  =  0 


(A22) 


apohu  +  ikhP  =  0 


(A23) 


cp08v  -l-  =  g 8p 


(A24) 
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After  eliminating  the  functions  8 u  and  8v,  we  have  two  equations  for  the  functions 
8p  and  8 P , 


o28p  4-  — 


(A26) 


(A27) 


Finally,  we  obtain  a  second  order  ordinary  differential  equation  for  8P, 


C2(r)-^f  +  +  C 0(2)8?  =  0, 


(A28) 
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(A30) 


(A31) 


Once  this  equation  can  been  solved,  all  the  other  three  functions  can  been  obtained 
from  8?.  8p  can  been  obtained  from  equation  (A27),  8 u  from  equation  (A23)  and 
8v  from  equation  (A24).  These  equations  hold  for  the  light  and  heavy  fluid 
separately.  They  are  coupled  together  by  the  boundary  conditions  at  their  interface; 
namely,  that  the  pressure  P  and  normal  velocity  are  continuous  to  first  order  in  8v. 
This  determines  the  growth  rate  a. 


To  illustrate  our  results,  we  assume  our  system  is  a  stiffened  polytropic  gas23 
The  specific  internal  energy  for  a  stiffened  polytropic  gas  is 
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(  P  +  yP, ) 

(  *Y  —  1  )p 


(A32) 


where  7  in  a  dimensionless  constant  and  Ps  is  a  constant  with  the  dimension  of  pres¬ 
sure.  The  unperturbed  gases  are  in  exponentially  stratified  isothermal  equilibrium: 

Pa  =  Pa0expOfl(z  -  2 intfc) )  (A33) 

Pi  =  PioexpOM?  -  2 intfc))  (A  34) 


where  ZjntfC  is  the  2  position  of  the  unperturbed  interface,  (3  =  and 

,1(P  +  P ,)  ,T 
c  =  (— —  1  ■ 

We  suppose  that  the  position  -q  of  the  perturbed  interface  is  given  by 
■>1  ~  ^intfc  =  A/m,exp(<xr)sin(Jbr)  ,  where  or  is  the  growth  rate  of  the  surface  and  k  is 
the  wave  number  of  the  perturbation  and  is  the  amplitude  of  the  perturbation  at 
interface  at  time  1  =  0.  We  assume  that  the  problem  is  periodic  in  x  with  reflecting 
boundaries  at  ;  =  2bdry  .  where  Zbdry  is  the  top  (bottom)  of  the  computational 
domain. 

The  linearized  Euler  equations  (A15)-(A18)  become 

tx8p  +  ikpoiu  +  Po~!p"  +  PPoBv  =  0  (A  35) 

crpo8/<  =  —  ikiP  (A  36) 


s-  .  diP  , 

apobv  +  -jp  =  gop 

(Tip  +  3(P0  +  P,) 8v  +  y(P0  +  P,)Vkhu  +  -^)  =  0  , 

and  equation  (A28)  becomes 

d2iP  _  yg  diP 
dz 2  c 2  &2 

Now  assume 

8 P  ~  exp  [(a  +  P)(z  -  zintfc)) 

8p  ~  exp  [(a  +  (3)(z  -  zintfc)j, 


(A37) 

(A38) 

(A39) 

(A40) 

(A41) 
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while 


8 u,  8v  ~  exp 


(A42) 


From  equation  (A39), 
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We  impose  the  boundary  conditions  8 v(z  =  Zbdry)  =  0  •  The  solution  of  the  linear¬ 
ized  Euler  equations  is  then 


8  P  =  - 


_ (*y  -  l)jg2  +  c2a 2 _ 

exp(a_ (Zjntfc  —  Zbdry))  CXp(tt  +  (Zjntfc  —  Zbdry)) 


PoAX 


exp(g-(z  -  Zbdry))  _  exp(a.(2  -  Zbdry)) 
(y  -  1)$  +  a-c 2  (y  -  l)g  +  a+c2 


exp(ar  +  ifcr)exp(3(z  -  Zjntfc)). 


(A44) 


8P  »  +  (7-1  1. 


8  u  = 
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8  P. 


(A45) 

(A46) 

(A47) 


Note  that  8v(z  =  ZintfC)  =  crAexp(af  +  ikx)  is  continuous  at  the  interface.  The  growth 
rate  a  is  determined  from  the  condition 


8 -  8 P.  = 


A  =  (p*  -  pa)gA 


(A48) 


which  expresses  the  continuity  of  the  pressure  at  the  interface 


P  +  (Zintfc  +  1|)  +  8P+  =  P -(Zjntfc  +  t|)  +  8 P.  .  (A49) 


Here  the  subscripts  +  (-)  denote  variables  evaluated  just  above  (below)  the  inter¬ 
face.  We  find 


(Pi  “  P a) 8  = 
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exp(aa_ (zintfc  “  Zbdry))  “  exp(afl+(zintfc  -  zbdry)) 
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FIGURE  CAPTIONS 

Fig.  4.1  The  location  in  x,  t  space  of  bubble  mergers.  ( H )  corresponds  to  the 
run  H,  and  ( L )  corresponds  to  run  L. 

Fig.  4.2H  bubble  radial  distribution  vs.  ln(r)  for  run  H  at  a  sequence  of  times. 

Fig.  4.2L  bubble  radial  distribution  vs.  ln(r)  for  run  L  at  a  sequence  of  times. 
Different  scales  are  used  between  the  ranges  from  0  to  60  and  from  60  to  5200  in  the 
y  direction  in  order  to  display  both  the  large  cluster  of  bubbles  near  r  =  ro  and  the 
remainder  of  the  distribution. 


Fig.  4.3  ln(number  of  bubbles)  vs.  time.  ( H )  corresponds  to  run  H,  and  (L) 
corresponds  to  run  L. 

Fig.  4.4  average  velocity  vs.  time.  ( H )  corresponds  to  run  H,  and  ( L ) 
corresponds  to  run  L.  The  derivative  of  the  velocity  versus  time  in  ( L )  is  about  7 
times  larger  than  the  derivative  in  (H) . 

Fig.  4.5  log  of  minimum ,  mean  and  maximum  bubble  radius  vs.  time.  ( H ) 

corresponds  to  run  H ,  and  (£.)  corresponds  to  run  L. 
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Fig.  4.6  ln(rmax/rm,a„)  vs.  time.  (H)  corresponds  to  run  H,  and  (L)  corresponds 
to  run  L. 


Fig.  4.7  minimum,  mean  and  maximum  bubble  height  vs.  time.  ( H )  corresponds 
to  run  H,  and  ( L )  corresponds  to  run  L. 

Fig.  4.8  radius,  height  cross  correlation  vs.  time.  ( H )  corresponds  to  run  H ,  and 
(L)  corresponds  to  run  L. 


Fig.  4.9  radius  and  height  correlation  vs.  time.  (H)  corresponds  to  run  H ,  and 
( L )  corresponds  to  run  L. 

Fig.  4.10  The  expected  scaled  radii  r~c  and  height  hc  of  bubbles  adjacent  to  large 
bubbles  are  plotted  as  a  function  of  time.  The  y  axis  is  in  units  of  dimensionless 
length.  This  plot  should  be  contrasted  to  Fig.  4.5  and  Fig.  4.7,  which  display  mean 
radius  and  height  without  restriction. 

Fig.  5.1  The  dependence  of  the  dimensionless  growth  rate  on  A  and  A/2.  The 
difference  between  two  adjacent  constant  contours  is  0.793. 

Fig.  5.2.  Plots  of  spike  velocity  and  bubble  velocity  versus  time  are  showm  with 
the  best  two  parameter  fit  to  equation  (21)  superimposed,  for  the  parameter  values 
D  =  2,  M2  =  0.5,  *y  =  1.4  .  The  numerical  results  are  obtained  by  using  a  80  by  640 
grid  in  a  computational  domain  1x8. 

Fig.  5.3  The  comparison  of  the  spike  velocity  and  the  spike  acceleration  of  the 
numerical  result  to  its  linear  and  large  time  asymptotic  behavior  for  D  =  2,  M2  =  .5 
and  y  =  1.4.  The  solid  lines  are  the  numerical  results  obtained  by  using  a  80  by  640 
grid  in  a  computational  domain  1x8. 

Fig.  5.4  The  plots  of  the  renormalization  factor  for  gravity  versus  the  Atwood 
ratio  for  ( a )  spike  and  ( b )  bubble.  For  comparison,  the  straight  lines  which  are  pro¬ 
portional  to  the  Atwood  ratio  are  also  plotted.  The  data  for  a  one  component, 
incompressible  fluid  (A  =  1,  Af2  =  0)  is  taken  from  reference  20. 
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Fig.  5.5  The  plot  of  the  spike  drag  coefficient  Cq  versus  the  Atwood  ratio  for 
A/2  =  .5  (  +  )  and  2  (X)  and  y  =  1.4.  Note  that  for  A  =  1,  Co  -  0  on  theoretical 
grounds  and  this  suggests  an  empirical  relation  Co  =  const  (1  —  A)  with  the  constant 
dependent  on  Af . 
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Fig.  5.6  Plots  of  Mach  number  (the  ratio  of  the  terminal  velocity  to  the  sound 
speed  of  the  heavy  gas)  of  (a)  spike  and  ( b )  bubble  versus  the  Atwood  ratio.  The 
high  density  ratio  systems  have  not  reached  the  terminal  velocity  region  at  the  end 
the  computation.  Therefore  the  terminal  velocity  can  not  be  determined  in  these  sys¬ 
tems.  The  spike  and  the  bubble  velocities  at  the  last  computation  step  are  supersonic 
(relative  to  the  sound  speed  of  the  heavy  gas)  in  the  D  =  10  and  100  cases.  In  the 
incompressible  case  (M2  =  0),  |vtenD/c|  =  0  by  definition. 

Fig.  5.7  A  plot  of  the  constant  c\  versus  the  square  root  of  the  Atwood  ratio  for 
D  =  2  and  3,  M2  =  .5  (  +  ),  2  (X),  and  y  =  1.4.  The  values  of  c\  calculated  from  the 
two  (*  and  A)  dimensional  incompressible  ( M 2  =  0,  A  =  1)  theories  are  also  shown. 
The  value  for  (*)  comes  from  reference  21  and  (A)  comes  from  reference  22. 

Fig.  5.8  Plots  of  the  interface  position,  density  and  pressure  contours  for 

D  =  1.5,  M2  =  .5,  a  =  .015,  y  =  1.4  in  a  computation  domain  1x6  with  a  40  by 

240  grid.  Only  the  upper  two  thirds  of  the  computational  region  is  shown  in  the  plot 
because  nothing  of  interest  occurs  in  the  remainder  of  the  computation,  (a)  The 
interface  position  for  successive  time  steps.  ( b )  The  density  contour  plot,  (c)  The 
pressure  contour  plot. 

Fig.  5.9  Plots  of  the  interface  position,  density  and  pressure  contours  for 
D  =  10,  A/2  =  .5,  a  =  .015,  y  =  1.4  in  a  computation  domain  1X4  with  a  40  by  160 
grid,  (a)  The  interface  position  for  successive  time  steps.  ( b )  The  density  contour 
plot,  (c)  The  pressure  contour  plot. 

Fig.  5.10  Plots  of  the  interface  position,  density  and  pressure  contours  for 

D  =  100,  M2  =  .5,  a  =  .015,  y  =  1.4  in  a  computation  domain  lx  10  with  a  20  by 

200  grid.  Only  the  upper  four  fifths  of  the  computational  region  is  shown  in  the  plot 
because  nothing  of  interest  occurs  in  the  remainder  of  the  computation,  (a)  The 
interface  position  for  successive  time  steps.  ( b )  The  density  contour  plot,  (c)  The 
pressure  contour  plot. 

Fig.  5.11  Plots  of  the  interface  position,  density  and  pressure  contours  for 
D  =  2,  M2  =  2,  a  =  .015,  y  =  1.4  in  a  computation  domain  lxg  with  a  80  by  640 
grid.  Only  the  middle  portion  of  the  computational  region  is  shown  in  the  plot 
because  nothing  of  interest  occurs  in  the  remainder  of  the  computation. 


Fig.  5.12  Plots  of  the  interface  position,  density  and  pressure  contours  for 
D  =  100,  M2  =  2,  a  =  .015,  y  =  1.4  in  a  computation  domain  1X6  with  a  20  by  120 
grid,  (a)  The  interface  position  for  successive  time  steps,  (b)  The  density  contour 
plot,  (c)  The  pressure  contour  plot. 

Fig.  5.13  Plots  of  the  interface  position,  density  and  pressure  contours  of  a  bub¬ 
ble  merger  for  D  =  2,  M2  =  2,  y  =  1.4  in  a  computation  domain  1x4  with  a  75  by 
300  grid,  (a)  The  interface  position  for  successive  time  steps.  ( b )  The  density  con¬ 
tour  plot,  (c)  The  pressure  contour  plot. 

Fig.  5.14  Plots  of  the  interface  position,  density  and  pressure  contours  of  a  bub¬ 
ble  merger  for  D  =  10,  M2  =  2,  *y  =  1.4  in  a  computation  domain  1x4  with  a  75  by 
300  grid,  (a)  The  interface  position  for  successive  time  steps,  (h)  The  density  con¬ 
tour  plot,  (c)  The  pressure  contour  plot. 

Fig.  6.1  Positions  of  the  front  at  a  sequence  of  time  steps  for  D  =  10,  M2  -  0.5 
and  a  =  0.015  by  using  40  by  160  and  10  by  40  grids.  The  finer  grids  resolve  the 
more  detailed  structure  of  the  interface. 

Fig.  6.2  A  comparison  of  the  position  of  the  spike  and  the  bubble  for  D  ~  10, 
M2  =  0.5  and  a  =  0.015  using  (A)  10  by  40,  (B)  20  by  80,  (C)  40  by  160  and  (H)  80 

by  320  grids.  Plot  a  is  the  position  of  the  spike.  Plot  b  is  the  position  of  the  bubble. 

\ 

For  comparison,  the  linear  theory  is  also  shown.  As  can  be  seen,  most  of  the  compu¬ 
tation  is  out  of  the  linear  regime. 

Fig.  6.3  A  comparison  of  the  spike  and  the  bubble  velocity  and  position  for  the 
linear  theory  and  the  full  Euler  equations.  Here  D  =  10,  M 2  =  0.5  and  a  -  0.001.  A 
40  by  160  grid  was  used  in  a  computational  region  of  size  1  by  4. 

Fig.  6.4  A  comparison  of  different  initial  starting  amplitudes  for  D  -  10, 
M2  =  0.5,  a  =  0.004  and  A  -  0.01  using  a  40  by  320  grid  in  a  computational  region 
of  size  1  by  4.  The  time  plot  is  terminated  at  t  =  4  because  for  a  =  0.004,  the  one 
way  reflected  waves  in  the  heavy  gas  arrive  at  the  interface  shortly  after  that  and  the 
agreement  between  these  two  runs  does  not  and  should  not  persist  after  this  time. 

Fig.  6.5  A  comparison  of  different  redistribution  frequencies  for  D  -  10, 
M 2  =  0.5  and  a  =  0.01  with  a  40  by  160  grid  in  a  computational  region  of  size  1  by 
4.  /  =  1,  15,  100  corresponds  to  redistributing  the  front  every  1,  15,  100  time  steps. 
Plot  a  is  for  the  position  of  the  spike.  Plot  b  is  for  the  interface.  Compare  to  frame 


Fig  6.6  A  comparison  of  boundary  effects  for  a  domain  of  shape  1X4  (A)  and 
1x8  (B)  with  10  by  40  and  10  by  80  grids  respectively.  Here  D  =  10,  M 2  =  0.5  and 
a  =  0.015.  The  one  way  reflected  wave  corresponds  to  a  signal  from  the  nearest 
boundary.  The  two  way  reflected  wave  corresponds  to  a  signal  starting  from  the 
interface.  Plot  a  is  for  the  velocity  of  the  bubble.  Plot  b  is  for  the  position  of  the 
bubble.  For  distinguishing  the  boundary  effects  from  the  grid  effects,  a  40  by  160 
grid  curve  (A’)  in  the  computation  domain  1X4  is  also  plotted  as  a  dotted  line  in  the 
plot  a.  The  wiggles  which  are  common  to  A  and  B  (and  missing  in  A')  are  grid 
effects,  while  the  divergence  between  A  and  B,  marked  by  the  arrival  of  the  one  way 
and  two  way  reflected  signals  are  pure  boundary  effects. 


Fig.  4.2L  bubble  radial  distribution  vs.  ln(r)  for  run  L  at  a  sequence  of  times.  Different 
scales  are  used  between  the  ranges  from  0  to  60  and  from  60  to  5200  in  the  y  direction 
in  order  to  display  both  the  large  cluster  of  bubbles  near  r  =  r0  and  the  remainder  of 
the  distribution. 
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Fig.  4.7  minimum,  mean  and  maximum  bubble  height  vs.  time.  (H)  corresponds  to  run  H, 
and  ( L )  corresponds  to  run  L. 
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Fig.  4.10  The  expected  scaled  radii  ~c  and  height  hc  of  bubbles  adjacent  to  large  bubbles 
are  plotted  as  a  function  of  time.  Hie  y  axis  is  in  units  of  dimensionless  length.  This 
plot  should  be  contrasted  to  Fig.  4.5  and  Fig.  4.7,  which  display  mean  radius  and  height 
without  restriction. 
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Fig.  5.2  Plots  of  spike  velocity  and  bubble  velocity  versus  time  are  shown  with  the  best 
two  parameter  fit  to  equation  (21)  superimposed,  for  the  parameter  values 
D  =  2,  M2  =  0.5,  y  =  1.4  .  The  numerical  results  are  obtained  by  using  a  80  by  640 
grid  in  a  computational  domain  1x8. 
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Fig.  S.3  The  comparison  of  the  spike  velocity  and  the  spike  acceleration  of  the  numeri¬ 
cal  result  to  its  linear  and  large  time  asymptotic  behavior  for  D  -  2,  M7  -  .5  and 
•y  «  1.4.  The  solid  lines  are  the  numerical  results  obtained  by  using  a  80  by  640  grid  in 
a  computational  domain  1x8. 


Fig.  5.4  \The  plots  of  th$  renormalization  factor  for  gravity  versus  the  Atwood  ratio  for 
(a)  spike  and  (b)  bubble.  For  comparison,  the  straight  lines  which  are  proportional  to 
the  Atwood  ratio  are  also  plotted.  The  data  for  a  one  component,  incompressible  fluid 
(A  *=  1,  A/2  =  0)  is  taken  from  reference  20. 


Fig.  5.5  The  plot  of  the  spike  drag  coefficient  CD  versus  the  Atwood  ratio  for 
M 2  «  .5  (+)  and  2  (x)  and  *y  =  1.4.  Note  that  for  A  «  1,  Cfl  =  0  on  theoretical 
grounds  and  this  suggests  an  empirical  relation  Q>  =  const  (1  -  A)  with  the  constant 
dependent  on  M. 
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Fig.  5.6  Plots  of  Mach  number  (the  ratio  of  the  terminal  velocity  to  the  sound  speed  of 
the  heavy  gas)  of  (a)  spike  and  (b)  bubble  versus  the  Atwood  ratio.  The  high  density 
ratio  systems  have  not  reached  the  terminal  velocity  region  at  the  end  the  computation. 
Therefore  the  terminal  velocity  can  not  be  determined  in  these  systems.  The  spike  and 
the  bubble  velocities  at  the  last  computation  step  are  supersonic  (relative  to  the  sound 
speed  of  the  heavy  gas)  in  the  D  =  10  and  100  cases.  In  the  incompressible  case 
(A/1  =  0),  |vtenn  /c|  =  0  by  definition. 
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Fig.  5.7  A  plot  of  the  constant  versus  the  square  root  of  the  Atwood  ratio  for  D  =  2 
and  3,  M 2  *  .5  (+),  2  (x),  and  y  =  1.4.  The  values  of  cj  calculated  from  the  two  (• 
and  A)  dimensional  incompressible  (M2  =  0,  A  =  1)  theories  are  also  shown.  The  value 
for  (*)  comes  from  reference  21  and  (A)  comes  from  reference  22. 
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Fig.  5.9  Plots  of  the  interface  position,  density  and  pressure  contours  for 
D  =  10,  M1  =  .5,  a  =  .015,  y  =  1.4  in  a  computation  domain  1x4  with  a  40  by  160 
grid,  (a)  The  interface  position  for  successive  time  steps.  (6)  The  density  contour  plot, 
(c)  The  pressure  contour  plot. 


Fig.  5.10  Plots  of  the  interface  position,  density  and  pressure  contours  for 
D  =  100,  Ml  -  .5,  a  -  .015,  y  =  1.4  in  a  computation  domain  1x10  with  a  20  by  200 
grid.  Only  the  upper  four  fifths  of  the  computational  region  is  shown  in  the  plot 
because  nothing  of  interest  occurs  in  the  remainder  of  the  computation,  (a)  The  inter* 
face  position  for  successive  time  steps.  ( b )  The  density  contour  plot,  (c)  The  pressure 
contour  plot. 
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Fig.  5  .11  Plots  of  the  interface  position  for  succuss  time  steps  for 
D  -  2,  Af:  *  2,  j4  =  .015,  y  =  1.4  in  a  computation  domain  1x8  with  80  by  640  grids. 
Only  the  middle  part  the  computational  region  is  shown  because  nothing  of  interest 
occurs  in  the  remainder  of  the  computation. 
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Fig.  5.12  Rots  of  the  interface  position,  density  and  pressure  contours  for 
D  -  100,  M2  =  2,  a  =  .015,  *y  =  1.4  in  a  computation  domain  1x6  with  a  20  by  120 
grid,  (a)  The  interface  position  for  successive  time  steps.  ( b )  The  density  contour  plot, 
(c)  The  pressure  contour  plot. 
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merger  for  D  -  2,  M2  =  2,  7  =  1.4  in  a  computation  domain  1x4  with  a  75  by  300 
grid,  (a)  The  interface  position  for  successive  time  steps.  ( b )  The  density  contour  plot, 
(c)  The  pressure  contour  plot. 


Fig.  5.14  Plots  of  the  interface  position,  density  and  pressure  contours  of  a  bubble 
merger  for  D  =  10,  M2  ~  2,  -y  *  1.4  in  a  computation  domain  1x4  with  a  75  by  300 
grid,  (a)  The  interface  position  for  successive  time  steps.  ( b )  The  density  contour  plot, 
(c)  The  pressure  contour  plot. 
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The  fronts  for  40  by  160  grids: 


The  fronts  for  10  by  40  grids: 


Fig.  6.1  Positions  of  the  front  at  a  sequence  of  time  steps  for  D  =  10,  M 2  -  0.5  and 
a  -  0.015  by  using  40  by  160  and  10  by  40  grids.  The  finer  grids  resolve  the  more 
detailed  structure  of  the  interface. 


Zspike 


Fig.  6.2  A  comparison  of  the  position  of  the  spike  and  the  bubble  for  D  *  10, 
M2  «  0.5  and  a  =  0.015  using  (A)  10  by  40,  (B)  20  by  80,  (C)  40  by  160  and  (E)  80  by 
320  grids.  Plot  a  is  the  position  of  the  spike.  Plot  b  is  the  position  of  the  bubble.  For 
comparison,  the  linear  theory  is  also  shown.  As  can  be  seen,  most  of  die  computation 
is  out  of  the  linear  regime. 
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Fig.  6.3  A  comparison  of  the  spike  and  the  bubble  velocity  and  position  for  die  linear 
theory  and  the  full  Euler  equations.  Here  D  -  10,  M :  =  0.5  and  a  «  0.001.  A  40  by 
160  grid  was  used  in  a  computational  region  of  sue  1  by  4. 


Fig.  6.4  A  comparison  of  different  initial  starting  amplitudes  for  D  -  10,  M 2  -  0.5, 
a  ■  0.004  and  A  =  0.01  using  a  40  by  320  grid  in  a  computational  region  of  size  1  by  4. 
The  time  plot  is  terminated  at  t  *=  4  because  for  a  —  0.004,  the  one  way  reflected  waves 
in  the  heavy  gas  arrive  at  the  interface  shortly  after  that  and  the  agreement  between 
these  two  runs  does  not  and  should  not  persist  after  this  time. 


-69- 


Vbubble 


Fig.  6.6  A  comparison  of  boundary  effects  for  a  domain  of  shape  1x4  (A)  and  1x8  (B) 
with  10  by  40  and  10  by  80  grids  respectively.  Here  D  -  10,  M 2  =  0.S  and  a  =  0.015. 
The  one  way  reflected  wave  corresponds  to  a  signal  from  the  nearest  boundary.  The 
two  way  reflected  wave  corresponds  to  a  signal  starting  from  the  interface.  Plot  a  is  for 
the  velocity  of  the  bubble.  Plot  b  is  for  the  position  of  the  bubble.  For  distinguishing 
die  boundary  effects  from  the  grid  effects,  a  40  by  160  grid  curve  (A')  in  the  computa¬ 
tion  domain  1x4  is  also  plotted  as  a  dotted  line  in  the  plot  a.  The  wiggles  which  are 
common  to  A  and  B  (and  missing  in  A’)  are  grid  effects,  while  the  divergence  between 
A  and  B,  marked  by  the  arrival  of  the  one  way  and  two  way  reflected  signals  are  pure 
boundary  effects. 
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